home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.7 / thin / thin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  25.9 KB  |  868 lines
  1. /* 
  2.  * thin.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* THIN:        program performs thinning on image 
  10.  *              usage: thin inimg outimg [-k K] [-n MAXITER] [-d DISPLAY] [-I] [-L]
  11.  *              The <display> argument shows result of each iteration
  12.  *              of thinning in "line-printer" style.
  13.  *                              
  14.  */
  15.  
  16. #include <stdio.h>
  17. #include <stdlib.h>
  18. #include <math.h>
  19. #include <string.h>
  20. #include <images.h>
  21. #include <tiffimage.h>          /* tiff file format info */
  22. extern void print_sos_lic ();
  23.  
  24. #define DFLTMAXK 3              /* dflt maximum thinning kernel size */
  25. #define MAXMAXK 21              /* max of the maximum thinning kernel size */
  26. #define MAXITER 20              /* maximum number of iterations */
  27. #define MIN0RUN 5               /* minimum run of zeros to run-length code */
  28. #define MAXDISPLAY 40           /* maximum length of line for display */
  29.  
  30. unsigned char OFF, ERASED, ON, PON;  /* values of pixels:
  31.                                       * * - ERASED value is for non-anchor
  32.                                       * * - ERASED  + 1 is for anchor
  33.                                       * * - ERASED increments +2 each iteration */
  34.  
  35. unsigned char **image;          /* input/output image */
  36. struct point imgSize;           /* image size */
  37. long ySizeM1;                   /* y length minus 1 */
  38. long **xRun;                    /* no., then x locns of 1/0 runs for each y */
  39.  
  40. long peel0 (unsigned char *, unsigned char *, long, long *, long *, long *);
  41. long peel (unsigned char *, unsigned char *, long, long *, long *);
  42. long ksize (long, long, long);
  43. long sqron (long, long, long);
  44. long getring (long, long, long, unsigned char *);
  45. long thinring (unsigned char *, long, int *);
  46. long chkconnect (unsigned char *, unsigned char *, long);
  47. long anchor (unsigned char *, long);
  48. long erasesqr (long, long, long, int, long *);
  49. long display (long);
  50. long usage (short);
  51. long input (int, char **, long *, long *, long *, long *);
  52.  
  53. main (argc, argv)
  54.      int argc;
  55.      char *argv[];
  56. {
  57.   register long x, y,           /* image coordinates */
  58.     i, k;                       /* sidelength of current thinning kernel */
  59.   Image *imgIO;                 /* structure for I/O images */
  60.   long maxK,                    /* max. sidelength of thinning kernel */
  61.     change[MAXMAXK],            /* no. erasures for each mask size */
  62.     nMaxIter,                   /* maximum number of iterations */
  63.     nIter,                      /* no. iterations */
  64.     nChange,                    /* no. thinning operations on iteration */
  65.     displayFlag,                /* display results after each iter if =1 */
  66.     invertFlag,                 /* invert input image before processing */
  67.     nONs,                       /* total ONs in original image */
  68.     nErased;                    /* cumulative no. ERASED in image */
  69.   unsigned char *ring;          /* ring of pixels on perimeter of kxk sqr */
  70.   unsigned char *side;          /* array of 8-connected side pixels */
  71.  
  72.   if ((input (argc, argv, &maxK, &nMaxIter, &displayFlag, &invertFlag)) < 0)
  73.     return (-1);
  74.  
  75.   imgIO = ImageIn (argv[1]);
  76.   image = imgIO->img;
  77.   imgSize.x = ImageGetWidth (imgIO);
  78.   imgSize.y = ImageGetHeight (imgIO);
  79.   if (imgSize.y > MAXDISPLAY)
  80.     displayFlag = 0;
  81.   ySizeM1 = imgSize.y - 1;
  82.  
  83. /* invert image */
  84.   if (invertFlag) {
  85.     for (y = 0; y < imgSize.y; y++)
  86.       for (x = 0; x < imgSize.x; x++)
  87.       image[y][x] = 255 - image[y][x];
  88.   }
  89.  
  90.   if ((xRun = (long **) calloc (imgSize.y, sizeof (long))) == NULL) {
  91.     printf ("CALLOC: xRun -- not enough memory -- sorry\n");
  92.     return (-1);
  93.   }
  94.  
  95.   if ((ring = (unsigned char *) malloc (4 * (maxK - 1))) == NULL) {
  96.     printf ("not enough memory -- sorry\n");
  97.     return (-2);
  98.   }
  99.   if ((side = (unsigned char *) malloc (maxK)) == NULL) {
  100.     printf ("not enough memory -- sorry\n");
  101.     return (-3);
  102.   }
  103.  
  104. /* initialize ON/OFF/ERASED values */
  105.   OFF = 0;
  106.   ERASED = 1;
  107.   ON = 255;
  108.   PON = ON - 1;
  109.  
  110. /* zero image borders */
  111.   for (y = 0; y < imgSize.y; y++)
  112.     image[y][0] = image[y][imgSize.x - 1] = OFF;
  113.   for (x = 0; x < imgSize.x; x++)
  114.     image[0][x] = image[ySizeM1][x] = OFF;
  115.  
  116.   for (k = 0; k < MAXMAXK; k++)
  117.     change[k] = 0;
  118.  
  119. /* on first iteration, perform thinning and accumulate information
  120.  * on x-runs */
  121.  
  122.   nChange = peel0 (ring, side, maxK, change, &nONs, &nErased);
  123.  
  124.   printf ("iteration  1:\t");
  125.   for (i = 3; i <= maxK; i++) {
  126.     printf (" %d) %3d;  ", i, change[i]);
  127.     change[i] = 0;
  128.   }
  129.   printf ("\n");
  130.  
  131.   if (displayFlag == 1)
  132.     display (nChange);
  133.   ERASED += 2;
  134.  
  135. /* iteratively convolve through image until thinned */
  136.  
  137.   for (nIter = 1, nChange = 1; (nChange > 0) && nIter <= nMaxIter; nIter++) {
  138.     nChange = peel (ring, side, maxK, change, &nErased);
  139.     printf ("iteration %2d:\t", nIter + 1);
  140.     for (i = 3; i <= maxK; i++) {
  141.       printf (" %d) %3d;  ", i, change[i]);
  142.       change[i] = 0;
  143.     }
  144.     printf ("\n");
  145.     if (displayFlag == 1)
  146.       display (nChange);
  147.     ERASED += 2;
  148.   }
  149.   printf ("%d thinned out of %d original ONs (%d%%)\n",
  150.           nONs - nErased, nONs, ((nONs - nErased) * 100) / nONs);
  151.   if (nIter >= nMaxIter && nChange != 0)
  152.     printf ("Nuts -- maximum iterations reached\n");
  153.  
  154.   for (y = 1; y < ySizeM1; y++)
  155.     for (x = 1; x < imgSize.x - 1; x++)
  156.       image[y][x] = (image[y][x] < PON) ? OFF : ON;
  157.  
  158. /* un-invert image */
  159.   if (invertFlag) {
  160.     for (y = 0; y < imgSize.y; y++)
  161.       for (x = 0; x < imgSize.x; x++)
  162.       image[y][x] = 255 - image[y][x];
  163.   }
  164.  
  165.   ImageOut (argv[2], imgIO);
  166.   return (0);
  167. }
  168.  
  169.  
  170.  
  171. /* PEEL0:       function performs first thinning iteration where information
  172.  *            on x-runs of ONs is accumulated
  173.  *                    usage: nChange = peel0 (ring, side, maxK, change, 
  174.  *                                                      &nOns, &nErased)
  175.  */
  176.  
  177. long
  178. peel0 (ring, side, maxK, change, nONs, nErased)
  179.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  180.      unsigned char *side;       /* array of 8-connected side pixels */
  181.      long maxK,                 /* max sidelength of thinning kernel */
  182.       *change,                  /* no. erasures for each mask size */
  183.       *nONs,                    /* total original ON pixels */
  184.       *nErased;                 /* cumulative no. ERASED in image */
  185. {
  186.   register long x, y,           /* image coordinates */
  187.     iXRun,                      /* index of runs in x */
  188.     xSizeM1,                    /* imgSize.x minus 1 */
  189.     k;                          /* sidelength of current thinning kernel */
  190.   long onRun,                   /* flag = 1 for run of 1s; 0 for 0s */
  191.     permOn,                     /* permanent ON flag set to 1 if so */
  192.     kM1,                        /* k minus 1 */
  193.     nChange;                    /* no. thinning operations on iteration */
  194.  
  195.   *nONs = *nErased = nChange = 0;
  196.   xSizeM1 = imgSize.x - 1;
  197.   for (y = 1; y < ySizeM1; y++) {
  198.     xRun[y] = (long *) calloc (imgSize.x + 1, sizeof (long));
  199.     xRun[y][0] = -MIN0RUN;
  200.     for (x = 1, iXRun = 1, onRun = 0; x < xSizeM1; x++) {
  201.       permOn = 0;
  202.       if (image[y][x] < ERASED) {
  203.         if (onRun == 1) {
  204.           onRun = 0;
  205.           xRun[y][iXRun++] = x - 1;
  206.         }
  207.         k = 0;
  208.       }
  209.       else {
  210.         if (onRun == 0) {
  211.           onRun = 1;
  212.           if ((x - xRun[y][iXRun - 1]) < MIN0RUN)
  213.             --iXRun;
  214.           else
  215.             xRun[y][iXRun++] = x;
  216.         }
  217.         (*nONs)++;
  218.         k = ksize (x, y, maxK);
  219.       }
  220.       kM1 = (k > 3) ? k - 1 : 3;
  221.       while (k >= kM1) {
  222.         if (sqron (x, y, k) == 0)
  223.           break;
  224.         if (getring (x, y, k, ring) == 1)
  225.           break;
  226.         if (thinring (ring, k, (int *) &permOn) == 1) {
  227.           if (chkconnect (ring, side, k) == 1) {
  228.             nChange++;
  229.             (change[k])++;
  230.             erasesqr (x, y, k, anchor (ring, k), nErased);
  231.             break;
  232.           }
  233.         }
  234.         --k;
  235.       }
  236.       if (permOn == 1)
  237.         image[y][x] = PON;
  238.     }
  239.     --iXRun;
  240.     if (iXRun % 2 != 0)
  241.       xRun[y][++iXRun] = x;
  242.     xRun[y][0] = iXRun;
  243.     xRun[y] = (long *) realloc (xRun[y], (sizeof (long)) * (iXRun + 1));
  244.   }
  245.   return (nChange);
  246. }
  247.  
  248.  
  249. /* PEEL:        function performs an iteration of thinning
  250.  *                usage: nChange = peel (ring, side, maxK, change, &nErased)
  251.  */
  252.  
  253. long
  254. peel (ring, side, maxK, change, nErased)
  255.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  256.      unsigned char *side;       /* array of 8-connected side pixels */
  257.      long maxK,                 /* max sidelength of thinning kernel */
  258.       *change,                  /* no. erasures for each mask size */
  259.       *nErased;                 /* cumulative no. ERASED in image */
  260.  
  261. {
  262.   register long x, y,           /* image coordinates */
  263.     xEnd, iXRun,                /* index of runs in x */
  264.     k;                          /* sidelength of current thinning kernel */
  265.   long noOnsInRow,              /* flag = 1 if no ONs in y-row (0 if ONs) */
  266.     permOn,                     /* permanent ON flag set to 1 if so */
  267.     kM1,                        /* k minus 1 */
  268.     nChange;                    /* no. thinning operations on iteration */
  269.  
  270.   nChange = 0;
  271.   for (y = 1; y < ySizeM1; y++) {
  272.     for (iXRun = 1, noOnsInRow = 1; iXRun <= xRun[y][0]; iXRun += 2) {
  273.       xEnd = xRun[y][iXRun + 1];
  274.       for (x = xRun[y][iXRun]; x <= xEnd; x++) {
  275.         permOn = 0;
  276.         if (image[y][x] < ERASED)
  277.           k = 0;
  278.         else
  279.           k = ksize (x, y, maxK);
  280.         kM1 = (k > 3) ? k - 1 : 3;
  281.         while (k >= kM1) {
  282.           if (sqron (x, y, k) == 0)
  283.             break;
  284.           noOnsInRow = 0;
  285.           if (getring (x, y, k, ring) == 1)
  286.             break;
  287.           if (thinring (ring, k, (int *) &permOn) == 1) {
  288.             if (chkconnect (ring, side, k) == 1) {
  289.               nChange++;
  290.               (change[k])++;
  291.               erasesqr (x, y, k, anchor (ring, k), nErased);
  292.               break;
  293.             }
  294.           }
  295.           --k;
  296.         }
  297.         if (permOn == 1)
  298.           image[y][x] = PON;
  299.       }
  300.     }
  301.     if (noOnsInRow == 1)
  302.       xRun[y][0] = -xRun[y][0];
  303.   }
  304.   return (nChange);
  305. }
  306.  
  307.  
  308.  
  309. /* KSIZE:       function determines k, where kxk is largest square
  310.  *            around (x,y) which contains all ON or ERASED
  311.  *                      usage: k = ksize (x, y, maxK)
  312.  */
  313.  
  314. long
  315. ksize (x, y, maxK)
  316.      long x, y,                 /* image coordinates */
  317.        maxK;                    /* maximum k value */
  318. {
  319.   register long xMask, yMask,   /* x,y mask coordinates */
  320.     xEnd, yEnd,                 /* end coord.s of square */
  321.     k;                          /* mask size */
  322.   long upHalf, downHalf,        /* half of mask below and above center */
  323.     xStart,                     /* x- start and end of square */
  324.     yStart;                     /* y- start and end of square */
  325.  
  326.   for (k = 4; k <= maxK; k++) {
  327.     if (k % 2 == 1)
  328.       downHalf = upHalf = (k - 3) / 2;
  329.     else {
  330.       upHalf = (k - 2) / 2;
  331.       downHalf = (k - 4) / 2;
  332.     }
  333.     xStart = x - downHalf;
  334.     xEnd = x + upHalf;
  335.     yStart = y - downHalf;
  336.     yEnd = y + upHalf;
  337.     for (yMask = yStart; yMask <= yEnd; yMask++)
  338.       for (xMask = xStart; xMask <= xEnd; xMask++)
  339.         if (image[yMask][xMask] < ERASED)
  340.           return (k - 1);
  341.   }
  342.   return (maxK);
  343. }
  344.  
  345.  
  346.  
  347.  
  348. /* SQRON:       function tests pixels in kxk square are already erased or 
  349.  *            for 3x3 if they can never be erased
  350.  *                      usage: flag = sqron (x, y, k)
  351.  *              flag = 0 if cannot erase any pixels in square
  352.  *                   = 1 if at least one pixel is still ON
  353.  */
  354.  
  355. long
  356. sqron (x, y, k)
  357.      register long x, y;        /* image coordinages */
  358.      long k;                    /* square sidelength of ring */
  359. {
  360.   register long xEnd, yEnd,     /* upper bounds of center erase area */
  361.     erasedP1;                   /* ERASED + 1 */
  362.   long upHalf, downHalf,        /* half of mask below and above center */
  363.     yStart, xStart;             /* bounds of center erase area */
  364.  
  365. /* check for 3x3 */
  366.   if (k == 3) {
  367.     if (image[y][x] == ON)
  368.       return (1);
  369.     else
  370.       return (0);
  371.   }
  372.  
  373. /* check center square */
  374.   if (k % 2 == 1)
  375.     downHalf = upHalf = (k - 3) / 2;
  376.   else {
  377.     upHalf = (k - 2) / 2;
  378.     downHalf = (k - 4) / 2;
  379.   }
  380.   xStart = x - downHalf;
  381.   xEnd = x + upHalf;
  382.   yStart = y - downHalf;
  383.   yEnd = y + upHalf;
  384.   erasedP1 = ERASED + 1;
  385.   for (y = yStart; y <= yEnd; y++)
  386.     for (x = xStart; x <= xEnd; x++)
  387.       if (image[y][x] >= erasedP1)
  388.         return (1);
  389.   return (0);
  390. }
  391.  
  392.  
  393.  
  394. /* GETRING:     function gets ring of pixels on perimeter of k-size square
  395.  *                usage: allOnes = getring (x, y, k, ring)
  396.  *              If ring includes only 1-value pixels, function returns 1,
  397.  *              otherwise 0.
  398.  */
  399.  
  400. long
  401. getring (x, y, k, ring)
  402.      register long x, y;        /* image coordinages */
  403.      long k;                    /* square sidelength of ring */
  404.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  405. {
  406.   register xEnd, yEnd,          /* x,y ends of square */
  407.     allOnes,                    /* =1 if all ones in ring; 0 otherwise */
  408.     i;
  409.   long upHalf, downHalf,        /* half of mask below and above center */
  410.     xStart,                     /* x- start and end of square */
  411.     yStart;                     /* y- start and end of square */
  412.  
  413.   if (k % 2 == 1)
  414.     downHalf = upHalf = (k - 1) / 2;
  415.   else {
  416.     upHalf = k / 2;
  417.     downHalf = (k - 2) / 2;
  418.   }
  419.   xStart = x - downHalf;
  420.   xEnd = x + upHalf;
  421.   yStart = y - downHalf;
  422.   yEnd = y + upHalf;
  423.  
  424.   allOnes = 1;
  425.   i = 0;
  426.   for (x = xStart, y = yStart; x <= xEnd; x++) {
  427.     if (image[y][x] < ERASED)
  428.       allOnes = 0;
  429.     ring[i++] = image[y][x];
  430.   }
  431.   for (y = yStart + 1, x = xEnd; y <= yEnd; y++) {
  432.     if (image[y][x] < ERASED)
  433.       allOnes = 0;
  434.     ring[i++] = image[y][x];
  435.   }
  436.   for (x = xEnd - 1, y = yEnd; x >= xStart; --x) {
  437.     if (image[y][x] < ERASED)
  438.       allOnes = 0;
  439.     ring[i++] = image[y][x];
  440.   }
  441.   for (y = yEnd - 1, x = xStart; y > yStart; --y) {
  442.     if (image[y][x] < ERASED)
  443.       allOnes = 0;
  444.     ring[i++] = image[y][x];
  445.   }
  446.  
  447. /* if this square is already at border, cannot go to larger square */
  448.   if (xStart <= 0 || yStart <= 0
  449.       || xEnd >= (imgSize.x - 1) || yEnd >= (ySizeM1))
  450.     return (0);
  451.  
  452.   return (allOnes);
  453. }
  454.  
  455.  
  456.  
  457. /* THINRING:    function makes decision to thin or not based on CNUM
  458.  *            and FNUM in perimeter ring
  459.  *                      usage: flag = thinring (ring, k, &permOn)
  460.  *              Flag = 1 if thinning conditions met, 0 otherwise.
  461.  */
  462.  
  463. long
  464. thinring (ring, k, permOn)
  465.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  466.      register long k;           /* square sidelength of ring */
  467.      int *permOn;               /* permanent ON flag set to 1 if so */
  468. {
  469.   register nRing,               /* no. pixels in ring */
  470.     cNum,                       /* connectivity number */
  471.     n, i;
  472.   long fNum,                    /* no. 1s on ring */
  473.     m;
  474.   long nOff,                    /* current run of OFF in ring */
  475.     phi0,                       /* maximum run of OFF in ring */
  476.     nFirstRun;                  /* no  OFF from ring[0] */
  477.   long lower, upper;            /* adjacent ring elements for cNum calc */
  478.  
  479.   nRing = 4 * k - 4;
  480.  
  481. /* calculate FNUM */
  482.   for (i = 0, fNum = 0; i < nRing; i++)
  483.     if (ring[i] >= ERASED)
  484.       fNum++;
  485.  
  486. /* calculate 4-connected run of 0s */
  487.   nOff = (ring[0] < ERASED) ? 1 : 0;
  488.   nFirstRun = (nOff == 1) ? 0 : -1;
  489.   for (i = 1, phi0 = 0; i < nRing; i++) {
  490.     if (ring[i] < ERASED)
  491.       nOff++;
  492.     else {
  493.       if (nOff > 0) {
  494.         phi0 = (nOff > phi0) ? nOff : phi0;
  495.         if (nFirstRun == 0)
  496.           nFirstRun = nOff;
  497.         nOff = 0;
  498.       }
  499.     }
  500.   }
  501.   if (nOff > 0) {
  502.     if (nFirstRun > 0)
  503.       nOff += nFirstRun;
  504.     phi0 = (nOff > phi0) ? nOff : phi0;
  505.   }
  506.  
  507. /* CNUM */
  508. /* CNUM skipping corners */
  509.   for (i = 2, cNum = 0; i < nRing; i++) {
  510.     lower = (long) ring[i - 1];
  511.     if ((i % (k - 1)) == 0)
  512.       i++;                      /* skip the corner pixels */
  513.     upper = (long) ring[i];
  514.     if (upper >= ERASED && lower < ERASED)
  515.       cNum++;
  516.   }
  517.   if (ring[1] >= ERASED && ring[nRing - 1] < ERASED)
  518.     cNum++;
  519. /* CNUM at corners */
  520.   for (n = 1; n < 4; n++) {
  521.     m = n * (k - 1);
  522.     if (ring[m] >= ERASED)
  523.       if (ring[m - 1] < ERASED && ring[m + 1] < ERASED)
  524.         cNum++;
  525.   }
  526.   if (ring[0] >= ERASED && ring[1] < ERASED && ring[nRing - 1] < ERASED)
  527.     cNum++;
  528.  
  529. /* to thin or not to thin */
  530.   if (cNum == 1)
  531.     if (phi0 > (k - 2) && fNum > (k - 2))
  532.       return (1);
  533.  
  534. /* for 3x3, set flag for perm. ON pixel if connection, end, or cross pt */
  535. /* note that 2nd conditional for cross added 17-Mar-89 to make equivalent
  536.  * to THINWX; this leaves non-8-conn. jct if cross-jct with one branch
  537.  * just one pixel, but it's consistent */
  538.   if (k == 3)
  539.     if (cNum > 1 || fNum <= 1 || (cNum == 0 && fNum == 4))
  540.       *permOn = 1;
  541.  
  542.   return (0);
  543. }
  544.  
  545.  
  546.  
  547. /* CHKCONNECT:  function checks connectivity to see if current erasures
  548.  *            combined with past erasures will destroy connectivity
  549.  *                      usage: flag = chkconnect (ring, side, k)
  550.  *              Function returns flag = 0 if connectivity destroyed,
  551.  *              flag = 1 if connectivity retained.
  552.  */
  553.  
  554. /* if corner or its nbr is ON, then that is value of side, otherwise OFF */
  555. #define CORNER(ISIDE,ICORNER,ICORNERSIDE,ICORNEROTHER) \
  556.         if (ring[ICORNER] >= PON || ring[ICORNEROTHER] >= ERASED) \
  557.               side[ISIDE] = ON; \
  558.         else side[ISIDE] = OFF
  559.  
  560. /* for a NW corner, check if corner is anchor point, and if so set to ON */
  561. #define CORNERNW(ISIDE,ICORNER,ICORNERSIDE,ICORNEROTHER) \
  562.         if (ring[ICORNER] >= ERASED + 1 \
  563.             || ring[ICORNEROTHER] >= ERASED) \
  564.               side[ISIDE] = ON; \
  565.         else side[ISIDE] = OFF
  566.  
  567. long
  568. chkconnect (ring, side, k)
  569.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  570.      unsigned char *side;       /* array of 8-connected side pixels */
  571.      long k;                    /* square sidelength of ring */
  572. {
  573.   register long i, nSide,       /* no. pixels along 8-connected side */
  574.     nSideM1,                    /* nSide minus 1 */
  575.     iSide,                      /* side array index */
  576.     nRing;                      /* no. pixels in ring */
  577.   long anON,                    /* an ON run along side */
  578.     ONandERASE,                 /* run of ON and n ERASEs */
  579.     N1;                         /* no. ones in nbrhood */
  580.  
  581.   nRing = 4 * k - 4;
  582.   nSide = k;
  583.   nSideM1 = nSide - 1;
  584.  
  585. /* calculate N1 */
  586.   for (i = 0, N1 = 0; i < nRing; i++)
  587.     if (ring[i] >= PON)
  588.       N1++;
  589.   if (N1 == 0)
  590.     return (0);
  591.  
  592. /* check connectivity of west side */
  593.   i = 3 * (k - 1);
  594.   CORNER (0, i, i + 1, i - 1);
  595.   for (i = i + 1, iSide = 1; i < nRing; i++, iSide++)
  596.     side[iSide] = ring[i];
  597.   CORNERNW (nSideM1, 0, nRing - 1, 1);
  598.  
  599.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  600.     if (side[i] >= PON) {
  601.       if (ONandERASE == 1)
  602.         return (0);
  603.       else
  604.         anON = 1;
  605.     }
  606.     else if ((side[i] == ERASED || side[i] == ERASED + 1)
  607.              && anON == 1)
  608.       ONandERASE = 1;
  609.     else if (side[i] < ERASED)
  610.       anON = ONandERASE = 0;    /* off */
  611.   }
  612.  
  613. /* check connectivity of north side */
  614.   i = k - 1;
  615.   CORNER (0, i, i - 1, i + 1);
  616.   for (i = i - 1, iSide = 1; i >= 0; --i, iSide++)
  617.     side[iSide] = ring[i];
  618.   CORNERNW (nSideM1, 0, 1, nRing - 1);
  619.  
  620.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  621.     if (side[i] >= PON) {
  622.       if (ONandERASE == 1)
  623.         return (0);
  624.       else
  625.         anON = 1;
  626.     }
  627.     else if ((side[i] == ERASED || side[i] == ERASED + 1)
  628.              && anON == 1)
  629.       ONandERASE = 1;
  630.     else if (side[i] < ERASED)
  631.       anON = ONandERASE = 0;    /* off */
  632.   }
  633.  
  634. /* check connectivity of east side */
  635.   i = k - 1;
  636.   CORNER (0, i, i + 1, i - 1);
  637.   for (i = i + 1, iSide = 1; iSide < nSideM1; i++, iSide++)
  638.     side[iSide] = ring[i];
  639.   CORNER (nSideM1, 2 * (k - 1), 2 * (k - 1) - 1, 2 * (k - 1) + 1);
  640.  
  641.   for (i = 0, anON = 0, ONandERASE = 0; i < nSide; i++) {
  642.     if (side[i] >= PON) {
  643.       if (ONandERASE == 1)
  644.         return (0);
  645.       else
  646.         anON = 1;
  647.     }
  648.     else if ((side[i] == ERASED || side[i] == ERASED + 1)
  649.              && anON == 1)
  650.       ONandERASE = 1;
  651.     else if (side[i] < ERASED)
  652.       anON = ONandERASE = 0;    /* off */
  653.   }
  654.   return (1);
  655. }
  656.  
  657.  
  658.  
  659. /* ANCHOR:      function returns 1 if core is exposed on NW side
  660.  *                    usage: anchor (ring, k)
  661.  */
  662.  
  663. long
  664. anchor (ring, k)
  665.      unsigned char *ring;       /* ring of pixels on perimeter of kxk sqr */
  666.      register long k;           /* square sidelength of ring */
  667. {
  668.   register long nRing,          /* no. pixels in ring */
  669.     i;
  670.  
  671.   nRing = 4 * k - 4;
  672.   for (i = 0; i < k; i++)
  673.     if (ring[i] >= ERASED)
  674.       return (0);
  675.   for (i = 3 * (k - 1) + 1; i < nRing; i++)
  676.     if (ring[i] >= ERASED)
  677.       return (0);
  678.   return (1);
  679. }
  680.  
  681.  
  682.  
  683. /* ERASESQR:    function erases square contained within square perimeter
  684.  *                    usage: erasesqr (x, y, k, anchor, &nErased)
  685.  *              If the core is an anchor, then the pixels are erased
  686.  *              to ERASED + 1, otherwise they are erased to ERASED.
  687.  *              For kxk > 3x3 erasure, PON pixels (permanent ON for 3x3)
  688.  *              are erased; and if an anchor point (ERASED + 1) can be
  689.  *              erased to a non-anchor point by a larger k (I don't know
  690.  *              if this can actually happen... I'm pretty sure it cannot
  691.  *              and have found on a fingerprint it does not, but have not
  692.  *              thought it out yet.) then it is erased to ERASED.
  693.  */
  694.    /* want to erase if ERASED +1, not just ON */
  695.  
  696. long
  697. erasesqr (x, y, k, anchor, nErased)
  698.      register long x, y;        /* image coordinages */
  699.      long k;                    /* square sidelength of ring */
  700.      register int anchor;       /* 1 if core is NW pt of NW-SE diag; else 0 */
  701.      long *nErased;             /* no. of erased */
  702. {
  703.   register long xEnd, yEnd;     /* upper bounds of center erase area */
  704.   long upHalf, downHalf,        /* half of mask below and above center */
  705.     yStart,                     /* bounds of center erase area */
  706.     xStart;
  707.  
  708. /* erase for 3x3 */
  709.   if (k == 3) {
  710.     if (image[y][x] == ON) {
  711.       (*nErased)++;
  712.       image[y][x] = (anchor == 0) ? ERASED : ERASED + 1;
  713.     }
  714.   }
  715. /* erase for kxk > 3x3 */
  716.   else {
  717.     if (k % 2 == 1)
  718.       downHalf = upHalf = (k - 3) / 2;
  719.     else {
  720.       upHalf = (k - 2) / 2;
  721.       downHalf = (k - 4) / 2;
  722.     }
  723.     xStart = x - downHalf;
  724.     xEnd = x + upHalf;
  725.     yStart = y - downHalf;
  726.     yEnd = y + upHalf;
  727.     for (y = yStart; y <= yEnd; y++) {
  728.       for (x = xStart; x <= xEnd; x++) {
  729.         if (image[y][x] >= ERASED + 1) {
  730.           if (image[y][x] >= PON)
  731.             (*nErased)++;
  732.           image[y][x] = (anchor == 0) ? ERASED : ERASED + 1;
  733.         }
  734.       }
  735.     }
  736.   }
  737.   return (0);
  738. }
  739.  
  740.  
  741.  
  742. /* DISPLAY:     function displays results after each iteration of thinning
  743.  *                    usage: display (nChange)
  744.  */
  745.  
  746. long
  747. display (nChange)
  748.      long nChange;              /* no. thinning operations on iteration */
  749. {
  750.   long y, x;                    /* image coordinates */
  751.   char c;
  752.  
  753.   for (y = 0; y < imgSize.y; y++) {
  754.     for (x = 0; x < imgSize.x; x++) {
  755.       if (image[y][x] == ON)
  756.         printf ("* ");
  757.       else if (image[y][x] == PON)
  758.         printf ("O ");
  759.       else if (image[y][x] == ERASED)
  760.         printf ("e ");
  761.       else if (image[y][x] == ERASED + 1)
  762.         printf ("a ");
  763.       else
  764.         printf ("  ");
  765.     }
  766.     printf ("\n");
  767.   }
  768.   if (nChange > 0) {
  769.     printf ("Enter <CR> to continue next iteration.\n");
  770.     scanf ("%c", &c);
  771.   }
  772.   return (0);
  773. }
  774.  
  775.  
  776.  
  777. /* USAGE:       function gives instructions on usage of program
  778.  *                    usage: usage (flag)
  779.  *              When flag is 1, the long message is given, 0 gives short.
  780.  */
  781.  
  782. long
  783. usage (flag)
  784.      short flag;                /* flag =1 for long message; =0 for short message */
  785. {
  786.  
  787. /* print short usage message or long */
  788.   printf ("USAGE: thin inimg outimg [-k K] [-n MAXITER] [-d DISPLAY] [-I] [-L]\n");
  789.   if (flag == 0)
  790.     return (-1);
  791.  
  792.   printf ("\nthin performs iterative thinning of binary objects\n");
  793.   printf ("in input image to produce skeleton image with values\n");
  794.   printf ("OFF (0) anf ON (255)\n\n");
  795.   printf ("ARGUMENTS:\n");
  796.   printf ("    inimg: input image filename (TIF)\n");
  797.   printf ("   outimg: output image filename (TIF)\n\n");
  798.   printf ("OPTIONS:\n");
  799.   printf ("     -k K: window size for kxk mask (k >= 3, default = %d)\n", DFLTMAXK);
  800.   printf (" -n NITER: maximum number of iterations (default max = %d)\n", MAXITER);
  801.   printf ("       -d: to display results of each iteration (< 40x40 image)\n");
  802.   printf ("       -I: invert input image before processing\n");
  803.   printf ("       -L: print Software License for this module\n");
  804.  
  805.   return (-1);
  806. }
  807.  
  808.  
  809. /* INPUT:       function reads input parameters
  810.  *                  usage: input (argc, argv, &maxK, &nMaxIter, &displayFlag)
  811.  */
  812.  
  813. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  814.  
  815. long
  816. input (argc, argv, maxK, nMaxIter, displayFlag, invertFlag)
  817.      int argc;
  818.      char *argv[];
  819.      long *maxK,                /* max. sidelength of thinning kernel */
  820.       *nMaxIter,                /* max. no. iterations */
  821.       *displayFlag,             /* display results after each iter if =1 */
  822.       *invertFlag;              /* invert input image before processing */
  823.  
  824. {
  825.   long n;
  826.   static char *kString = "-k";
  827.   static char *nString = "-n";
  828.   static char *dString = "-d";
  829.   static char *IString = "-I";
  830.   static char *dashString = "-";
  831.  
  832.   if (argc == 1)
  833.     USAGE_EXIT (1);
  834.   if (argc == 2)
  835.     USAGE_EXIT (0);
  836.  
  837.   *maxK = DFLTMAXK;
  838.   *nMaxIter = MAXITER;
  839.   *displayFlag = 0;
  840.   *invertFlag = 0;
  841.  
  842.   for (n = 3; n < argc; n++) {
  843.     if (strcmp (argv[n], kString) == 0) {
  844.       if (++n == argc || argv[n][0] == '-')
  845.         USAGE_EXIT (0);
  846.       *maxK = atol (argv[n]);
  847.     }
  848.     else if (strcmp (argv[n], nString) == 0) {
  849.       if (++n == argc || argv[n][0] == '-')
  850.         USAGE_EXIT (0);
  851.       *nMaxIter = atol (argv[n]);
  852.     }
  853.     else if (strcmp (argv[n], dString) == 0)
  854.       *displayFlag = 1;
  855.     else if (strcmp (argv[n], IString) == 0)
  856.       *invertFlag = 1;
  857.     else if (strcmp (argv[n], "-L") == 0) {
  858.       print_sos_lic ();
  859.       exit (0);
  860.     }
  861.     else
  862.       USAGE_EXIT (0);
  863.   }
  864.   if (*maxK < 3)
  865.     USAGE_EXIT (1);
  866.   return (0);
  867. }
  868.